##########################################################################################

library(data.table)
library(optparse)
library("dndscv")

##########################################################################################

option_list <- list(
    make_option(c("--info_file"), type = "character"),
    make_option(c("--input_im_file"), type = "character"),
    make_option(c("--input_igc_file"), type = "character"),
    make_option(c("--input_dgc_file"), type = "character"),
    make_option(c("--out_path"), type = "character")
)

if(1!=1){
    
    info_file <- "~/20220915_gastric_multiple/dna_combinePublic/public_ref/combine/MutationInfo.combine.addMolecularSubType.rmMIX.tsv"
    input_im_file <- "~/20220915_gastric_multiple/dna_combinePublic/maf_public/All_use.IM.maf"
    input_igc_file <- "~/20220915_gastric_multiple/dna_combinePublic/maf_public/All_use.IGC.maf"
    input_dgc_file <- "~/20220915_gastric_multiple/dna_combinePublic/maf_public/All_use.DGC.maf"
    out_path <- "~/20220915_gastric_multiple/dna_combinePublic/finalPlot/revise/smgs"

}

parseobj <- OptionParser(option_list=option_list, usage = "usage: Rscript %prog [options]")
opt <- parse_args(parseobj)
print(opt)

info_file <- opt$info_file
input_im_file <- opt$input_im_file
input_igc_file <- opt$input_igc_file
input_dgc_file <- opt$input_dgc_file
out_path <- opt$out_path

##########################################################################################

dat_mut_info <- data.frame(fread(info_file))
dat_mut_im <- data.frame(fread(input_im_file))
dat_mut_igc <- data.frame(fread(input_igc_file))
dat_mut_dgc <- data.frame(fread(input_dgc_file))

## 只看NJMU
#dat_mut_info <- subset(dat_mut_info , From=="NJMU")

##########################################################################################

dat_mut_im$Class <- "IM"
dat_mut_igc$Class <- "IGC"
dat_mut_dgc$Class <- "DGC"
dat_mut <- rbind( dat_mut_im , dat_mut_igc , dat_mut_dgc )
dat_mut <- subset(dat_mut , Tumor_Sample_Barcode %in% dat_mut_info$Tumor)
print(dim(dat_mut))

##########################################################################################
### https://htmlpreview.github.io/?http://github.com/im3sanger/dndscv/blob/master/vignettes/dNdScv.html
##   sampleID chr      pos ref mut
## 1 Sample_1   1   871244   G   C
## 2 Sample_1   1  6648841   C   G
## 3 Sample_1   1 17557072   G   A

result <- c()
for( classN in unique(dat_mut$Class) ){
    ## 提取突变类型
    dat_mut_use <- subset( dat_mut , Class == classN )
    ## 同一患者的相同突变合并
    mutations <- unique(dat_mut_use[,c("Tumor_Sample_Barcode" , "Chromosome" , "Start_position" , "Reference_Allele" , "Tumor_Seq_Allele2")])
    colnames(mutations) <- c("sampleID" , "chr" , "pos" , "ref" , "mut")
    ## 开算
    dndsout <- dndscv(mutations , max_muts_per_gene_per_sample = Inf, max_coding_muts_per_sample = Inf)
    sel_cv = dndsout$sel_cv
    sel_cv$Class <- classN
    result <- rbind( result , sel_cv )
}

out_name <- paste0( out_path , "/compute_dn_ds.csv" )
write.csv( result , out_name , row.names = F , quote = F )